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FOREWORD 


Physical  processes  often  cannot  be  accurately  modeled  as  a  purely  deterministic 
mathematical  process  because  of  the  existence  of  random  aspects  of  their  behavior.  This 
unpredictable  portion  is  often  statistically  modeled  as  Gaussian  white  noise.  However,  noise 
processes  are  not  always  characterized  by  the  properties  of  Gaussian  white  noise.  For  example,  a 
projectile  whose  motion  is  disturbed  by  air  or  water  turbulence  or  a  self-propelled  projectile  with 
unevenly  mixed  fuel  may  possess  long-term  slowly  decreasing  dependencies  that  may  not  be 
well  represented  by  white  noise.  This  research  concerns  an  optimal  estimator  for  parameters  and 
states  of  systems  driven  by  another  type  of  noise,  known  as  fractional  Gaussian  noise  (FGN) 
processes.  These  stochastic  processes  can  model  systems  containing  long-term,  slowly 
decreasing  time-correlated  random  disturbances. 

This  report  examines  an  estimator  for  an  unknown  parameter  in  a  model  represented  by  a 
form  of  FGN-driven  stochastic  differential  equation.  A  simulation  of  fractional  Brownian 
motion  (FBM)  process  and  a  state  model  driven  by  FGN  are  also  developed  for  the  purpose  of 
testing  the  estimator.  FBM  is  a  stochastic  process  that  corresponds  to  FGN.  This 
correspondence  is  discussed.  The  estimator  is  tested  using  the  simulations  of  FBM  and  the  state 
model  driven  by  FGN.  An  alternate  algorithm  for  estimating  the  unknown  parameter  is  derived 
in  this  report.  Also,  estimating  the  single  parameter  that  controls  the  time  correlation  of  FBM  is 
discussed. 

Time  and  funding  constraints  did  not  permit  further  study  to  compare  the  methodologies 
reported  here  with  more  conventional  methods  for  a  specific  application. 

The  research  involved  in  the  report  was  in  support  of  the  Navy  In-house  Laboratory 
Independent  Research  (IR)  Program  at  the  Naval  Surface  Warfare  Center,  Dahlgren  Division 
(NSWCDD).  This  report  has  been  reviewed  by  Dr.  Edward  J.  Wegman,  Dunn  Professor, 
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INTRODUCTION 


The  motion  of  a  projectile  or  torpedo  traveling  through  the  air  or  water  may  be  governed 
by  a  differential  equation  whose  variable,  commonly  called  its  state,  represents  speed.  Random 
disturbances  affecting  the  flight  or  trajectory  of  the  object  may  be  caused  by  air  or  water 
turbulence,  or,  in  the  case  of  a  self-propelled  projectile,  the  disturbance  could  have  originated 
from  an  uneven  fuel  mixture.  Many  physical  processes,  such  as  these,  often  cannot  be  accurately 
modeled  as  a  purely  deterministic  mathematical  process  because  of  the  existence  of  random 
aspects  of  their  behavior.  This  unpredictable  portion  of  the  phenomena,  often  known  as  noise,  is 
frequently  statistically  modeled  as  Gaussian  white  noise.  However,  noise  processes  are  not 
always  characterized  by  the  properties  of  white  noise.  Thus,  Gaussian  white  noise,  which  by 
definition  assumes  time  independence  may  not  foiTn  an  accurate  model.  Furthermore,  a  constant, 
such  as  one  representing  friction  divided  by  projectile  mass  may  represent  an  unknown 
parameter  in  the  difterential  equation  (see  Reference  1,  p.  154).  This  investigation  concerns 
optimal  estimators  for  parameters  and  states  of  systems  driven  by  a  type  of  noise  that  possesses 
long-term,  slowly  decreasing  time  dependencies.  This  noise  is  known  as  fractional  Gaussian 
noise  (FGN),  and  it  may  represent  systems  with  slowly  decreasing  dependencies,  such  as  the 
example  just  described,  more  accurately  than  traditional  white  noise  models. 

The  use  of  stochastic  processes  in  the  form  of  stochastic  differential  equations 
(differential  equations  disturbed  by  random  noise)  driven  by  white  noise  for  modeling  physical 
phenomena  has  been  extensively  studied  for  the  case  where  the  equations  were  driven  by 
Gaussian  white  noise.  Parametric  estimators  have  been  developed  for  parameters  in  some  forms 
of  stochastic  differential  equations  driven  by  this  type  of  noise.  This  report  examines  an 
estimator  that  was  described  in  Reference  2  that  can  be  used  to  estimate  an  unknown  parameter 
of  a  model  represented  by  a  form  of  FGN-driven  stochastic  differential  equation.  A  simulation 
of  fractional  Brownian  motion  (FBM)  process  and  a  state  model  driven  by  FGN  was  also 
developed  for  the  purpose  of  testing  the  estimator.  FBM  is  a  stochastic  process  that  corresponds 
to  FGN.  The  estimator  was  tested,  using  the  simulations  of  FBM  and  the  state  model  driven  by 
FGN.  An  alternate  estimator  was  also  derived  in  this  report. 

Another  parameter  that  is  very  crucial  in  any  models  involving  FBM  is  the  single 
parameter  known  as  H  that  characterizes  FBM  itself.  The  higher  the  value  of  H,  the  greater  the 
time  dependency  of  the  FBM.  In  this  research,  H  is  assumed  to  be  within  the  interval  (0.5,  1.0) 
although  it  can  be  defined  for  all  H  within  the  interval  (0.0, 1.0).  The  interpretation  for  processes 
where  H  ^  1  is  unclear  (Reference  3,  p.  943).  After  investigation,  an  algorithm  was  found  in 
literature  for  estimating  the  single  parameter  needed  to  characterized  FBM  and  its  corresj)onding 
Gaussian  noise  process.  This  parameter  determines  the  time  correlation  of  the  FBM  and  FGN. 
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BACKGROUND 


This  section  provides  some  basic  definitions  and  mathematical  results  needed  for 
developing  an  optimal  estimator  for  parameters  and  states  of  FGN-driven  systems.  These 
definitions  and  results  can  also  help  in  understanding  the  reasons  for  potential  applications  of  this 
work. 


FBM  DEFINITION 

The  most  important  background  item  is  the  definition  of  FBM,  because  this  process  is  the 
originator  of  all  processes  in  this  research. 

DEFINITION:  Let  B  =  {B(t):  t  e  R}  be  a  standard  Brownian  motion  (BM)  process  and  R  the  set 
of  real  numbers.  Fractional  Brownian  motion  (FBM),  symbolized  by  Bh,  given  H  e  (0.5,  1)  is 
defined  as  follows: 


"  rT^TTTT^  JC* "  )dB(s)  + 

1  +  0.5  ) 

/  —  oo 

!  H-0  5  (1) 

J  It-  sr  °-^dB(s)) 

0 

where  T  is  the  gamma  function,  and  s  varies  from  negative  infinity  to  t. 

Hence,  FBM  is  a  stochastic  integral  of  a  detemiinistic  function  with  respect  to  standard 
Brownian  motion. 

For  H  =  0.5,  FBM  is  standard  BM.  Hence,  FBM  can  be  thought  of  as  a  generalization  of 
standard  BM.  FBM  is  a  zero  mean  normally  distributed  process  with  the  following  covariance 
function  (R): 


2 


(2) 


V  _-r(2-2H)cos(pH) 

"  ;tH(2H-l)  ’  ^ 


H 


(3) 


The  distribution  of  Bpi  also  has  the  property  of  self-similarity,  which  can  formally  be 
defined  for  FBM  as  follows: 


d  H 

{Bh  (at):t  e  R}  =  a  e  R };  d  =  distribution  (4) 
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Equation  (4)  defines  FBM  self-similarity  and  states  that  the  distribution  of  FBM  looks  the  same 
regardless  of  the  scale  factor  in  time.  The  self-similarity  property  is  also  commonly  known  as 
the  fractal  property,  and  therefore  the  distribution  of  Br  can  be  said  to  possess  the  fractal 
property.  One  can  see  that  by  setting  H  =  0.5,  BM  also  possesses  this  fractal  property.  This 
property  is  useful  because  some  physical  processes,  or  the  random  noise  part  of  the  physical 
processes,  appear  to  behave  as  fractals. 

An  example  of  a  physical  system  possessing  random  noise  that  behaves  like  FBM  is 
random  errors  in  communication  channels  as  found  in  Reference  3.  The  errors  in  these  channels 
may  appear  as  groups  of  bursts,  and  this  group  of  bursts  is  itself  grouped  in  bursts,  giving  the 
appearance  of  the  self-similarity  property. 

The  FBM  can  also  be  defined  by  Equation  (1)  for  H  e  (0,  0.5).  However,  for  modeling 
long-term  dependencies,  H  is  considered  to  be  between  0.5  and  1. 


FGN  DEHNITION 

FGN  (symbolized  by  Wr)  is  considered  to  be  the  derivative  of  FBM,  just  as  Gaussian 
white  noise  is  considered  to  be  the  derivative  of  standard  Brownian  motion.  In  this  sense,  FGN 
corresponds  to  FBM  and  vice  versa.  Yet,  analogous  to  the  fact  that  Brownian  motion  is  almost 
surely  (a.  s.)  not  differentiable,  we  have  the  following  theorem: 

THEOREM  1  (Reference  4):  Fractional  Brownian  motion  is  a.  s.  not  differentiable  with  respect 
tot. 

The  significance  of  this  theorem  is  to  show  that  in  the  strict  mathematical  sense,  FGN  (along 
with  Gaussian  white  noise)  does  not  exist  as  a  true  stochastic  process.  Yet,  it  is  still  desirable  to 
use  them  as  processes  because  of  the  useful  properties  they  possess:  FGN's  long-term  slowly 
decreasing  time  dependencies  and  1/f  properties  (described  in  the  next  paragraph)  and  Gaussian 
white  noise  being  time  independent.  These  properties  approximate  physical  systems. 
Furthermore,  physical  phenomena  often  behave  as  if  they  were  derivatives  of  FBM  or  Brownian 
motion  in  spite  of  the  fact  that  such  derivatives  do  not  truly  exist.  See  Reference  3  for  a 
description  of  how  to  define  FGN  in  terms  of  an  operator  instead  of  a  derivative  of  FBM. 

FGN  can  be  thought  of  as  a  stochastic  process  defined  to  be  normally  distributed  with 
zero  mean  and  the  covariance  function  (R),  which  is  defined  in  Equation  (5),  and  spectral  density 
function  (Swr)’  defined  in  Equation  (6): 


R(t,s)  =  R^  (x)=Vh  H  (2H-l)ltl 


2H-2 


H 


H 


V  -r(2-2H)cos(pH)  T:  =  t-s 
“  7tH(2H-l) 


^w  (w)=lwl^  ^^;w6R,W9i0 


(5) 


(6) 
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Equation  (5)  (the  covariance  function)  shows  that  FGN  exhibits  tlie  long-term  dependency 
property  of  having  slowly  decreasing  con-elations.  Because  -1.0  <  2H  -  2  <  0.0,  the  correlation  is 
high  for  small  t,  which  decreases  slowly  with  increasing  t  since  2H  -  2  >  -1.0.  Moreover,  since 
-1.0  <  1  -  2H  <  0.0,  the  spectral  density  function  (Equation  (6))  shows  that  FGN  has  low- 
frequency,  slowly  decreasing  spectral  power  for  the  same  reason  that  FGN  has  high  correlation 
for  small  t  (decreasing  slowly  with  increasing  x).  Processes  with  such  spectral  densities  are 
known  as  1/f  processes.  Along  with  being  normally  distributed,  these  unique  properties  give 
FGN  its  potential  for  applications  as  mathematical/statistical  models  of  physical  processes. 

Note  that  by  setting  x  to  zero.  Equation  (5)  shows  that  FBM  has  infinite  variance.  Also, 
notice  that  since  -1.0  <  1  -  2H  <  0, 


fSyr,  (w)dw  =  jlwl^  ^^dw  =  oo 


Figure  1  contains  graphs  showing  sample  paths  of  standard  Brownian  motion  and 
Gaussian  white  noise  (H  =  0.5).  Figure  2  shows  graphs  of  a  sample  path  of  FBM  and  FGN, 
where  H  =  0.85.  These  graphs  are  presented  merely  to  show  pictorially  examples  of  the  noise 
processes  being  studied.  They  are  not  part  of  the  findings  in  this  investigation. 


SPECTRAL  DENSITY 

Since  spectral  density  can  be  inteipreted  as  the  average  spectral  power  per  unit  frequency 
(Reference  5,  pp.  92-93),  Equation  (7)  shows  that  FGN  has  infinite  spectral  power.  Similarly, 
Gaussian  white  noise  has  infinite  variance  and  infinite  average  spectral  power.  No  natural 
process  can  possess  such  properties.  Hence,  in  the  physical  sense  as  well  as  in  the  mathematical 
sense,  neither  FGN  nor  Gaussian  white  noise  can  exist.  However,  both  of  these  pseudo¬ 
processes  can  form  approximations  of  many  real-life  phenomena  where  the  spectral  power 
density  is  non-zero  for  a  wide  range  of  frequencies. 

Gaussian  white  noise  is  well  known  in  the  literature  for  having  a  constant  spectral  density 
function.  Therefore,  it  is  often  useful  for  representing  processes  with  constant  spectral  power 
over  a  wide  range.  FGN  is  useful  for  representing  phenomena  in  an  approximate  manner  with 
low  frequency  power  that  is  high  near  zero  and  slowly  decreasing  over  a  wide  range  before 
becoming  zero.  Note  that  processes  with  the  long-term  slowly  decreasing  time  dependencies,  as 
given  by  Equation  (5)  (covariance)  used  for  "defining"  FGN,  can  be  shown  to  have  the  low 
frequency  spectral  power  as  also  given  in  Equation  (6). 

One  example  of  a  physical  process  that  behaves  as  a  FGN  is  discharge  from  a  river. 
These  discharges  may  tend  to  exhibit  clusters  of  high  and  low  periods  and  therefore  exhibit  long¬ 
term  dependencies.  As  a  second  example,  the  influence  of  sea  wave  action  on  the  accuracy  of 
projectiles  fired  from  weapons  on  ships  may  also  possess  the  long-term  time  dependency 
property  of  FGN.  Further  investigation  needs  to  be  conducted  regarding  this  second  example  to 
study  how  closely  it  can  be  approximated  by  FGN.  Taking  the  "derivative"  of  a  model  of  the 
errors  portion  of  the  communications  channel  (as  described  previously)  gives  a  third  example  of 
FGN.  Also,  the  effects  of  air  or  water  turbulence  on  projectiles  or  the  effects  of  unevenly  mixed 
fuel  could  be  candidates  for  FGN  representation.  This  was  described  in  the  introduction. 
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FIGURE  2.  RANDOM  NOISE  II  Oi  =  0.85) 
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Although  FBM  is  not  differentiable,  a  "pseudo-derivative"  can  be  approximated  by 
determining  the  ratio  of  an  increment  of  Bh  to  an  increment  of  time  and  thinking  of  it  as  a 
derivative  in  analyzing  mathematical  properties  of  the  model. 


MARTINGALE 

Another  important  concept  is  the  definition  of  a  martingale.  A  martingale  is  a  stochastic 
process  where  the  expectation  of  its  future  state  is  precisely  its  present  state.  The  precise 
definition  follows: 

DEFINITION:  Given  a  pair  (M(t),  a(t):  t  e  T)  where  M(t)  is  a  stochastic  process  that  is 
measurable  with  respect  to  a(t),  an  increasing  family  of  a-algebras,  then  the  pair  is  called  a 
martingale  if 


E  I  M  ( t )  I  <  CO  and 

(8) 

E[M(t)la(s)]  =  M(s)Vt>s 

The  most  common  example  of  such  a  proce.ss  is  standard  Brownian  motion  itself.  See 
Reference  2  p.  6  for  a  definition  of  a  a-algebra. 

A  martingale  has  two  useful  properties:  Future  expectations  can  be  described,  and  a 
martingale  lends  itself  to  well-defined  stochastic  integrals.  If  a  process  is  not  at  least  a 
generalization  of  martingales  known  as  a  local  martingale,  possible  problems  arise  in  defining 
integrals  of  another  stochastic  process  with  respect  to  this  original  non-martingale  process.  See 
Reference  5,  p.  234,  for  a  precise  definition  of  a  local  martingale. 


PROBLEM  STATEMENTS  AND  APPROACH 
OPTIMAL  ESTIMATOR  FOR  PARAMETER  0 


Developing  the  E.stimator 

The  first  problem  was  to  develop  an  optimal  estimator  for  the  parameter  0  in  the 
following  FGN-driven  stochastic  differential  equation. 

dX  (t)  =  0X(t)dt  +  W^j  (t)dt  (9) 

With  FGN  not  being  a  technically  true  stochastic  process,  perhaps  a  more  accurate  notation  in  the 
strict  mathematical  sense  for  the  same  equation  is  shown  in  Equation  (10): 

dX  (t)  =  0X(t)dt  +  dB„  (t)  (10) 
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Because  FBM  is  not  differentiable,  one  may  ask  what  is  the  meaning  of  this  stochastic 
differential  equation.  The  answer  is  that  Equation  (10)  is  a  symbolic  form  of  the  following 
integral  equation: 


X(t)  =  eJX(s)ds+BH(t)  (11) 

0 

Comparing  the  stochastic  differential  equation  with  the  random  noise  in  the  form  of  FGN  with 
the  equivalent  stochastic  integral  equation  with  its  random  noise  in  the  form  of  FBM  further 
reinforces  the  concept  that  FBM  and  FGN  correspond  to  one  another. 

The  stochastic  differential  equation,  Equation  (10),  with  an  unknown  parameter  6  is 
potentially  useful  in  modeling  physical  systems  whose  behavior  resembles  that  of  a  family  of 
differential  equations  indexed  by  the  parameter  6  and  whose  random  noise  behaves  like  FGN. 
Hence,  finding  an  estimate  of  the  unknown  parameter  9  is  one  method  of  selecting  which  of  the 
differential  equations  in  such  a  family  of  equations  best  approximates  the  system  being  modeled. 

Fitting  Equation  (10),  or  equivalently  Equation  (9),  into  the  projectile  example  given  in 
the  introduction,  6  represents  the  frictional  constant  divided  by  the  projectile  mass,  X  represents 
the  projectile  velocity,  and  Wh  represents  the  random  disturbance. 

Methods  have  been  established  to  derive  an  optimal  estimator  if  the  noise  process  is  a 
square- integrable  martingale.  However,  one  can  show  the  following  theorem: 

THEOREM  2  (Reference  2,  pp.  47-49):  Bh  is  not  a  local  martingale  for  H  e  (0.5, 1). 

Therefore,  Bh  cannot  be  a  martingale  for  H  e  (0.5, 1)  because  local  martingales  are  more  general 
than  martingales.  A  direct  proof  that  Bh  is  not  a  martingale  for  H  e  (0.5, 1)  was  also  presented  in 
Reference  2.  Hence,  a  method  of  estimating  the  parameter  in  Equation  (10)  needed  to  be 
developed.  As  already  noted,  when  H  =  0.5,  Bh  is  standard  Brownian  motion,  which  is  a 
martingale. 

One  of  the  methods  developed  in  Reference  2  for  estimating  0  was  based  upon  a  pseudo¬ 
least-squares  estimator  shown  in  Christopeit  (Reference  6).  In  Christopeit’s  "quasi-least-squares" 
method,  the  random  noise  is  still  assumed  to  be  a  martingale.  Otherwise,  however,  one  can  show 
that  the  model  represented  by  Equation  (10)  can  fit  into  the  model  given  in  Reference  6.  Note 
that  the  quasi-least- squares  method  was  derived  as  a  continuous  generalization  of  the 
conventional  (discrete)  least- squares  estimator.  A  conventional  least-squares  estimator  was 
derived  without  regard  to  random  noise  distribution.  Also,  the  only  places  where  the  martingale 
property  was  used  in  developing  the  quasi-least-squares  method  were  in  defining  stochastic 
integrals  and  in  detemining  asymptotic  properties  such  as  consistency.  Hence,  the  least-squares 
estimator  itself  will  be  mathematically  correct  as  long  as  the  terms  within  the  estimator, 
involving  stochastic  integrals,  are  well-defined. 
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The  modified  Christopeit's  “qua  si -least-squares”  solution  is  shown  in  Equation  (12): 


JX(s)dX  (s) 


9(t)=± 


Jx2(s)ds 


(12) 


See  Reference  2,  pp.  63-65,  for  a  detailed  explanation  of  this  formula. 

The  denominator  on  the  right-hand  side  of  Equation  (12)  is  merely  a  Riemann  or 
Lebesgue  integral.  However,  the  interpretation  of  the  integral  in  the  numerator  is  not  as  clear. 
The  form  of  the  model’s  stochastic  differential  equation  shows  that  the  dX  (s),  located  under  the 
integral  sign  in  the  numerator,  is  actually  equivalent  to  0X  (s)  ds  +  dBn  (s);  that  is. 


dX  (s)  =  0X(s)ds  +  dBj^  (s) 


(13) 


Hence,  the  integral  in  the  numerator  of  the  right-hand  side  implicitly  contains  two  integrals:  an 
integral  of  0X  (s)  with  respect  to  s,  which  can  be  inteipreted  as  a  Riemann  or  Lebesgue  integral, 
and  a  stochastic  integral  with  respect  to  Br,  which  is  in  the  following  foim 


jx(s)dBH(s) 


(14) 


These  facts  were  explained  in  Reference  2,  pp.  64-65.  Since  Br  is  neither  a  square-integrable 
martingale  or  a  local  martingale,  investigation  was  required  to  determine  whether  this  integral 
can  be  shown  to  be  well-defined.  An  argument  developed  in  Reference  2,  pp.  53-62,  shows  that 
this  integral  is  indeed  well  defined.  This  finding  was  essential  to  assure  the  estimator's 
numerator  is  meaningful. 

The  formula  for  the  modified  Christopeit's  quasi-least-squares  estimation  method  was 
coded  in  FORTRAN  for  running  on  a  VAX  computer  system. 


Testing  the  Estimator 

To  test  the  program,  a  method  for  simulating  the  X  given  parameter  6  must  be  found. 
The  algorithm  that  was  derived  in  this  research  was  in  the  form  of  a  Monte  Carlo  simulation  of 
the  solution  to  the  model  shown  in  Equation  (10).  The  solution  to  this  stochastic  differential 
equation  follows  for  X  (0)  =  0: 


X(t)  =  |e(‘-^^®dB^(s) 
0 


(15) 
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The  derivation  of  Equation  (15)  is  explained  in  Reference  2,  pp.  37-40.  Note  that  because  the 
integrand  of  this  integral  is  a  deterministic  function,  there  is  no  ambiguity  concerning  whether 
the  integral  is  well-defined. 

Because  X  is  an  integral  with  respect  to  FBM,  in  order  to  simulate  X  in  a  Monte  Carlo 
manner,  a  method  of  simulating  the  FBM  for  any  H  such  that  H  e  (0.5,  1)  must  also  be  found.  A 
Monte  Carlo  simulation  of  Bh  was  developed  as  a  multivariate  normally  distributed  random 
vector  of  dimension  1000  or  2000  where  each  component  Bh  (i)  (i  =  1,...,  n;  n  =  1000  or  2000) 
represented  a  different  time,  i  At,  of  Bh.  At  was  a  constant  to  be  read  in  as  an  input  representing 
the  desired  increment  in  time  between  successive  components  of  the  random  vector. 

Ideally,  one  desires  to  have  the  highest  possible  dimension  in  the  random  vector  because 
the  higher  the  dimension,  the  greater  the  amount  of  data,  which  improves  the  ability  to  test  for 
the  convergence  of  the  estimator.  However,  because  of  capacity  limitations  of  the  VAX 
computer,  either  1000  or  2000  data  points  are  near  the  maximum  amount  of  data  for  inputting 
into  the  computer  for  this  particular  project.  A  dimension  of  2000  was  used  for  analysis  in  single 
precision,  but,  because  double  precision  requires  more  memory,  a  dimension  of  only  1000  was 
used  for  double  precision  analysis.  The  mean  of  the  nonually  distributed  random  vector  was  set 
to  be  zero  because  FBM  by  definition  has  zero  mean.  The  covariance  matrix  was  designed  so 
that  each  entry  represented  a  covariance  of  the  FBM  between  two  of  the  times. 

Given  the  simulated  FBM  and  noting  that  X's  solution  was  an  integral  with  respect  to 
FBM,  X  was  simulated  as  a  numerical  approximation  of  the  integral  with  respect  to  the  simulated 
FBM.  A  numerical  approximation  of  a  Monte  Carlo  simulation  of  X  thus  was  obtained. 

Finally,  noting  that  Christopeit  used  the  martingale  property  for  determining  consistency 
and  knowing  that  Bh  lacks  this  property,  a  third  task  was  to  find  an  alternate  method  for 
determining  consistency.  After  extensive  study,  an  analytic  algorithm  has  not  yet  been  found. 
Therefore,  one  possible  method  may  be  to  empirically  test  the  consistency  using  the  simulated  X. 
The  empirical  testing  of  consistency  has  not  yet  been  performed.  Perfonning  an  empirical  test  of 
consistency  is  currently  difficult  because  the  VAX's  memory  capability  limits  the  dimension  of 
the  multivariate  random  vector  FBM  simulation  to  1000  or  2000  (as  previously  explained). 

An  alternate  estimator  is  also  derived  in  this  report. 


FBM  PARAMETER  H  ESTIMATION 

Another  problem  for  investigation  was  to  estimate  the  FBM  parameter  H,  assuming  that 
H  e  (0.5,  1).  In  the  Equation  (10)  model,  the  FBM,  Bh  (t),  was  assumed  to  be  unknown  for 
nonnegative  t  because  only  X  (t)  or  increments  of  X  (t)  were  assumed  to  be  observable. 
However,  in  some  applications,  Bh  (t),  may  possibly  lend  itself  to  observation  for  negative  t. 
Hence,  two  data  sources  may  be  used  to  estimate  H. 

The  most  accurate  way  is  to  estimate  H  directly  from  the  observed  data  for  Bh  (0  for 
negative  t.  If  this  data  is  unavailable,  the  alternate  method  is  to  first  estimate  0  and  then  solve  for 
the  Bh  (t)  in  the  integral  form  of  the  stochastic  differential  equation  model: 

B„(t)  =  X(t)-eJX(s)ds  (16) 

0 
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The  estimated  Bh  (t)  obtained  in  this  manner  can  then  be  used  to  estimate  H,  although  the 
accuracy  of  the  estimated  H  from  data  obtained  by  this  alternate  method  is  unknown  and  needs  to 
be  studied,  probably  by  empirical  methods. 

One  method  for  estimating  H  (given  the  FBM  data)  is  described  in  Reference  7.  This 
method  is  based  upon  a  least-squares  fit  to  the  incremental  Bh  data  given  by 

k 

j»k)  —  ^  [ABji  ((j+  i)At)  *“  m  (ABj^  (( j-i-  i)At))]^ 


where  j  is  any  positive  integer.  The  A  represents  the  increment  of  Bh  between  (j-fi)At  and 
(j-i-i-l)At.  H  is  estimated  as  being  one-half  (1/2)  the  slope  of  a  straight  line  fitted  using  the  least- 
squares  method  to  data  given  in  the  form  of  log  V  (j,  k)  versus  log  (k)  as  shown  in  Equation  (18). 

logV(j,k)=0.5Hlogk-t-C  (18) 


C  is  a  constant. 


KALMAN  FILTER  GENERALIZATION  FOR  X  ESTIMATOR 

Another  research  problem  also  dealt  with  the  same  FGN-driven  stochastic  differential 
^uation,  Equation  (9);  that  is,  the  same  equation  with  FGN  as  the  random  noise.  One  assumed 
in  this  problem  that  parameter  0  was  known  but  the  state  X  (t)  was  unknown.  Moreover,  an 
observation  model  with  the  random  component  being  the  traditional  Gaussian  white  noise  was 
assumed  to  be  available.  The  problem  was  to  investigate  the  possibility  of  generalizing  the 
Kalman  filter  into  an  estimator  that  can  estimate  the  state  X  in  spite  of  the  fact  that  X  was  driven 
by  FGN  instead  of  Gaussian  white  noise. 


ACHIEVEMENTS 


The  previous  section  defined  problems  that  must  be  considered  for  FGN  Estimation 
Theory.  Approaches  for  resolving  those  problems  were  also  examined.  This  section  details  the 
procedures  and  achievements. 


OPTIMAL  ESTIMATOR  FOR  PARAMETER  0 

As  stated  in  the  previous  section,  an  estimator  for  0  in  the  stochastic  differential  equation, 
Equation  (10)  was  developed  in  the  form  of  the  modified  Christo peit's  generalized  least-squares 
fit  (Equation  (12)).  The  coding  for  this  estimator  and  every  computer  program  in  this  research 
project  was  written  in  FORTRAN  using  the  VAX  computer  system.  This  estimator  was  checked 
by  using  the  simulation  of  X  that  was  a  function  of  the  Monte  Carlo  simulation  of  Bh. 
Reference  8  described  the  software  package  that  was  used  to  generate  a  normally  distributed 
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random  vector.  Theorem  2  showed  that  Bh  is  not  a  local  martingale  for  H  e  (0.5,  1).  Also, 
stochastic  integrals  involved  in  the  estimator  can  be  shown  to  be  well-defined  (see  previous 
section  or  Reference  2,  pp.  53-62,  pp.  64-65).  This  is  necessary  to  ensure  that  the  estimator  itself 
is  well-defined. 

Denoting  the  single  parameter  of  FBM  again  by  H  and  assuming  H  is  in  the  interval 
[0.5,1),  this  algorithm  worked  well  for  H  =  0.5  for  three  conditions: 

(1)  The  unknown  parameter  is  positive  and  less  than  or  equal  to  approximately  4.6. 

(2)  The  increment  of  time  is  set  to  0.001  units. 

(3)  The  array  is  allowed  to  hold  a  sequence  of  2000  different  times  (see  Table  1). 

The  numbers  as  shown  for  0  in  Table  1  were  chosen  merely  as  representations  of  possible 
parameters  and  are  not  of  any  particular  significance.  If  the  parameter  was  set  to  a  value  much 
larger  than  approximately  4.6,  the  program  aborted  because  X's  magnitude  was  too  large. 


TABLET.  ESTIMATION  OF  9 


VALUE  OF  H 

6=  1.0 

0  =  4.5 

0  =  -2.0 

0  »  4.6 

0.5 

e  =  0.09 

£  =  0.01 

e  =  0.08 

program  aborts 

0.75 

e  =  0.005 

e  =  0.1 

£=1.9 

It  was  noted  that  FBM  becomes  standard  Brownian  motion  when  H  =  0.5.  A  true  FBM 
that  is  not  merely  Brownian  motion  was  simulated  with  H  =  0.75.  In  this  case,  double  precision 
was  needed  for  the  program  to  work.  However,  double  precision  necessitated  greater  computer 
memoi7  than  single  precision.  Therefore,  the  array  had  to  be  decreased  back  to  1000  elements. 
Also,  the  increment  of  time  was  changed  to  0.01  units.  For  positive  valued  parameters  not 
exceeding  4.6,  reasonably  accurate  estimates  resulted.  However,  the  results  were  not  as 
favorable  for  the  negative  valued  parameter,  probably  beeause  for  the  negative  value  parameters, 
1000  different  time  elements  did  not  comprise  a  sufficiently  large  sampling  of  the  model. 
However,  in  double  precision,  the  VAX's  space  limitation  did  not  allow  an  array  much  larger 
than  1000.  The  FORTRAN  listing  of  the  double  precision  versions  of  the  simulator  and 
estimator  programs  is  shown  in  Appendix  A.  Table  1  also  gives  results  for  H  =  0.75.  Further 
investigation  of  the  handling  of  this  array  problem  is  needed  to  determine  alternative  possibilities 
for  testing  the  estimator. 

Another  method  for  simulating  the  FBM  for  use  in  obtaining  X  is  to  first  perform  a 
Monte  Carlo  simulation  of  independent  Gaussian  random  numbers,  which  can  be  used  to 
represent  dB  (s)  for  time  s.  Then  the  FBM  can  be  simulated  using  either  the  given  formula 
defining  Bh  (see  Equation  (1))  or  another  equation  for  FBM  restricted  to  time  within  a  closed 
and  bounded  set  [0,  T],  such  that  T  <  oo;  that  equation  is  shown  in  the  next  section.  Then  the  X 
can  be  simulated  in  the  same  manner  already  de.scnbed. 
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Refer  to  Equation  (12)  for  the  estimator.  Suppose  that  the  experimental  data  is  "close"  to 
being  continuous  and  analog,  enabling  readings  for  extremely  small  increments  of  time.  As 
shown,  the  estimator  involves  a  ratio  of  two  stochastic  integrals,  defined  to  be  the  quadratic  mean 
(q.  m.)  limits  of  summations  over  partitions  of  [0,  t].  This  convergence  is  analyzed  below. 

E  [X2  (s)]  is  a  finite  and  continuous  function  for  all  se[0,t]  for  finite  t  with  the  set  [0,  t] 
being  Lebesgue  measurable.  Therefore,  the  integral  in  the  denominator  can  be  treated  as  a 
Lebesgue  integral,  according  to  Reference  5,  p.  45,  instead  of  as  a  quadratic  mean  integral. 
Furthermore,  Equation  (15)  showed  X  as  an  integral  of  a  continuous  function  with  respect  to  a 
FBM,  which  is  the  sum  of  2  integrals  of  continuous  functions  with  respect  to  a  Brownian  motion. 
Brownian  motion  is  equivalent  to  a  process  that  is  a.  s.  sample  path  continuous  (Reference  9, 
pp.  66-67).  This  implies  that  X^  must  be  equivalent  to  an  a.  s.  continuous  process.  Thus,  the 
integral  can  further  be  assumed  to  be  a  Riemann  integral.  The  denominator  as  a  Riemann 
integral  is  analyzed  as  follows: 

Consider  the  partition  of  [0,  t],  {tg  =  0,  ti,...,tn  =  t}.  Let  I](t;6,n)  represent  the  numerical 
summation  approximating  the  integral  of  the  denominator.  Then  lli(t;q,n)  -  Ii(t;q,m)l->0  as 
and  max(At)— >0  according  to  Cauchy's  sequence  theorem.  This  numerical  integration 
can  be  repeatedly  performed  to  obtain  both  Ii(t;e,m)  and  Ii(t;e,n),  refining  the  partition  with 
each  iteration,  assuming  that  the  continuous  experimental  data  allows  such  refinement,  until 
lll(t;0,n)  -  Ii(t;6,m)l  no  longer  decreases. 

Since  the  numerator  is  not  equivalent  to  any  Lebesgue  or  Riemann  integral,  analyzing  the 
partition  of  [0,  t]  is  more  involved  than  analyzing  the  denominator.  Also,  an  estimated  value  for 
H  and  a  finite  range  of  possible  values  of  9  must  be  known.  For  this  case,  let  l2(t,9,n)  be  the 
function  approximating  the  integral  in  the  denominator.  Since  l2(t,9,n)  converges  to  the  actual 
integral  in  q.  m.  instead  of  in  the  ordinary  limit  (as  in  the  Lebesgue  or  Riemann  integral), 
Il2(t;9,n)-  l2(t;9,m)l-^0  cannot  be  claimed.  The  q.  m.  convergence  implies  that  VE[  Il2(t;0,n)- 
l2(t;9,m)|2]-^0.  However,  without  knowing  the  precise  value  for  9  and  given  only  a  single 
sample  path  of  data,  E[  Il2(t;9,n)-  l2(t;9,m)|2]  cannot  be  determined.  Thus,  the  values  for 
E[  Il2(t;9,n)-  l2(t;9,m)|2]  must  be  determined  for  a  known  range  of  possible  values  of  9  or  by 
substituting  the  already  estimated  value.  One  can  refine  the  partition  {to=0,  ti,...,  tn=t}  until 
VE[  Il2(t;9,n)-  l2(t;9,ra)|2]  no  longer  decreases  for  each  9  and  from  this  set  choose  the  one  with 
the  most  refined  partition.  Note  that  given  a  trial  9,  VE[  Il2(t;9,n)-  l2(t;9,m)|2]  can  be  obtained  by 
two  methods.  One  way  is  by  mathematical  derivation  using  Equation  (15)  for  X(t)  and  the 
properties  of  expectations.  The  other  is  by  Monte  Carlo  simulation  of  a  set  of  sample  paths  of 
l2(t;9,n)-  l2(t;9,m)  and  finding  the  sample  second  moment. 

The  puipose  of  this  description  of  convergence  is  only  to  describe  its  nature.  In  actuality, 
the  practical  way  to  check  the  partition  for  a  given  t  is  to  merely  select  a  partition  and  repetitively 
refine  it  until  the  estimated  9  stabilizes. 

An  alternate  estimator  and  its  derivation  are  discussed  below.  It  has  the  advantage  of  a.  s. 
convergence,  which  will  be  justified.  However,  its  disadvantage  is  that  the  value  of  the  H 
parameter  for  the  FBM  is  assumed  to  be  already  known,  while  the  estimator  just  described  does 
not  make  this  assumption. 
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FBM  PARAMETER  H  ESTIMATING  THE  STATE  X 

An  algorithm  from  Reference  7  that  is  an  estimator  for  the  value  of  the  H  parameter, 
given  FBM  data,  was  described  in  the  previous  section.  Future  plans  to  implement  this  estimator 
in  FORTRAN  and  to  test  it  using  the  simulated  FBM  have  been  proposed.  Also  proposed  are 
plans  to  test  this  estimator  for  its  performance  when  the  data  are  available  only  in  the  form  of  the 
simulated  X,  assuming  that  the  value  of  the  parameter  9  is  known. 


KALMAN  FILTER  GENERALIZATION  FOR  ESTIMATING  THE  STATE  X 

Finally,  considerable  research  has  been  performed  in  an  attempt  to  derive  a  generalization 
of  the  Kalman  filter.  That  filter  can  perform  an  optimal  estimate  of  the  state  X  for  the  model 
where  the  parameter  0  was  then  assumed  to  be  known  (see  Equation  (10)).  Recall  that  this 
problem  is  a  generalization  of  the  Kalman  filter  since  the  noise  is  FGN  and  not  Gaussian  white 
noise,  which  is  the  assumed  noise  of  a  state  model  to  be  estimated  by  the  Kalman  filter.  At  the 
present  time,  however,  such  a  generalized  filter  has  not  yet  been  worked  out.  Whether  such  a 
generalized  filter  can  be  mathematically  derived  cuiTently  remains  unknown. 


ALTERNATE  ESTIMATOR  FOR  PARAMETER  9 

Suppose  that  the  value  of  H  is  unknown  and  the  model  is  as  follows: 

dX  (t)  =  9X  (t)dt  +  dBjj  (t )  for  t>0  and 

dX(t)=dBH(t)  otherwise 

The  H  parameter  can  then  be  estimated  by  using  the  data  for  negative  time.  With  this  estimated 
H,  another  algorithm,  derived  in  Reference  2,  can  be  used  to  estimate  the  0.  This  algorithm  was 
shown  to  be  strongly  consistent.  In  other  words,  the  alternate  estimator  converges  a.  s.  to  the 
actual  0. 

Now  suppose  that  data  are  only  available  for  t  e  [0,T]  for  finite  T  and  that  the  value  of  H 
is  assumed  to  be  known.  The  following  alternate  algorithm  that  converges  a.  s.  as  finite  T 
becomes  infinitely  large  can  be  derived  for  the  model: 

dX  (t)  =  0X(t)dt  +dBj^  (t)  te[0,T] 


The  following  theorem  must  be  used  to  derive  this  new  estimator: 

THEOREM  3  (Reference  3,  pp.  951-952):  Given  the  closed  and  bounded  indexing  set  [0,t], 
there  exists  a  standard  Brownian  motion  process  Bx  for  te  [0,T]  such  that  for  the  FBM 
restricted  to  t  e  [0,T],  symbolized  by  Bhit.  the  following  transfoitnation  holds: 


Bx(t) 


t  S 

H~0.5 


r(i.5-H)  ^ 
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With  the  X  in  the  model,  let 


Y(t)  = 

This  implies  that: 

Y(t)  = 


r(1.5-H)  ^ 


f(s-ur-Hu°^-»dX(u) 


(22) 


1 


[0JsH-«-^d,  J(s-u)°-^-»  u«-^-»X(u)du  + 


r(  1.5 -H)  0 


Js  H  J(s  -  U  ^  U  0-5  H(}B  (u)] 

0  0 


(23) 


Theorem  3  implies  that  the  second  term  for  the  expression  of  Y  must  be  Brownian  motion. 
Hence, 


Y(t) 


1 


(1.5  -H) 


[6  Js»-o-5d  J(s-u)®-^““  X(u)  du+  (t)] 


(24) 


or 


Y(t)  =  0R(t)  +  B.j,(t) 


where 


(25) 


R(t) 


1 

r  ( 1.5-H) 


J3H-0.5 

0 


s 

dj  j(s  -  u ) u  X (u  )du 
0 


(26) 


The  formula  for  integration  by  parts  is  known  to  also  hold  true  for  integrals  of  a  stochastic 
process  with  respect  to  time  or  integrals  of  a  deterministic  function  with  respect  to  a  stochastic 
process  (see  Reference  9,  p.  88).  This  fact  is  needed  to  find  the  derivative  of  R;  R  is  needed  for 
using  the  estimator  that  is  now  being  derived.  Another  fact  used  to  find  dR/dt  is  Theorem  4: 

THEOREM  4:  Let  f  and  g  be  two  deterministic  and  differentiable  functions  in  [0,T].  Let  X  (not 
necessarily  the  same  X  of  the  model  in  this  report)  be  any  stochastic  process  with  a  covariance 
function  that  is  of  bounded  variation  in  the  set  [0,T]  x  [0,T].  Then 

t  t  t 

J  f (s)d[g(s)X(s)]  =J  f(s)X(s)g(s)ds  +  1  f (s)g(s)dX  (s)  (27) 

0  0  0 
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Proof:  By  the  stochastic  version  of  the  integration  by  parts  formula, 

f  t 

J  f(s)d[g(s)X(s)]  =  [f(s)g(s)X(s)]J)  -  Jg(s)X(s)Us)ds 


(28) 


Moreover,  by  using  the  stochastic  version  of  integration  by  parts  again, 

r  ‘  ‘ 

J  f(s)g(s)dX  (s)=  [f(s)g(s)X(s)]i  -  Jg(s)X(s)f(s)ds  -J  f(s)X(s)g(s)ds 


0 

Therefore, 

t 


(29) 


Jf(s)X(s)g(s)ds+ J  f(s)g(s)dX  (s)  =  [f(s)g(s)X(s)]o  -  J  g(s)X(s)f(s)ds 
^0  0 

(30) 

Note  that  both  the  left  and  right  sides  of  the  equation  in  the  theorem  statement  are  equal  to  the 
same  expression,  meaning  that  the  equation  must  be  true.  Therefore,  the  proof  of  Theorem  4  is 
now  complete. 

Also,  if  a  process,  say  Q,  is  either  sample  path  or  q.m.  differentiable. 


Jf(s)dQ  (s)=Jf(s)(^^-^)ds 

ds 


(31) 


where  the  derivative  is  either  q.m.  or  sample  path,  and  f  is  an  ordinary  deterministic  function. 
This  fact  can  be  proven  by  using  the  two  variations  of  the  stochastic  integral  by  parts  formula  to 
show  that  either  of  the  above  two  integrals  are  equal  to  the  same  thing. 

With  these  concepts.  Equation  (26)  for  R(t)  can  now  be  manipulated  to  give  a  form  that 
lends  itself  to  deriving  its  derivative.  Using  Leibniz  formula  for  derivatives  on  the  differential 
term, 


d  J(s-r)0-5-Hx(r)dr  = 

0  (32) 


S-8 

lim  [  1  (0.5  -  H  )(s  -  r  r^^-H  x(r  )dr  +  e  o^-h  (s_  g)  o.s-h  x  (s-e)]ds 

E->0  0  /  V  /J 


Note  that  the  above  integral  is  a  function  of  s  that  is  of  bounded  variation  and  thus  has  a 
derivative  almost  everywhere  (a.  e.)  in  s.  Hence,  the  traditional  sample  path  and  the  quadratic 
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mean  derivative  are  a.  e.  the  same  in  s.  In  analyzing  the  above  differential  term,  one  can  think  of 
it  as  a  sample  path  differential  that  obeys  the  conventional  laws  of  differential  calculus  a.  e. 
With  the  conventional  laws  of  calculus  and  uniform  convergence  in  mind  while  considering  ds, 
the  limit  of  d/ds  and  d/ds  of  the  limit  can  be  shown  to  be  equal  for  the  case  being  considered 
here.  Explanations  can  be  found  in  books  on  advanced  calculus  or  real  analysis  such  as 
Reference  10.  Making  the  substitutions  for  integration  by  parts  on  the  term  with  the  integral  in 
the  above  expression, 

u  =  X(r  )r0-5-H  du  =(  0.5  -  H  )X (r  )r  “®-^^dr  +r®-^"”dX(r)  (33) 

dv  =  (s- r)~°-^"“dr=»  v=  "q '  (s  -  r)  (34) 

Equation  (33)  for  du  is  merely  a  set  of  symbols  from  formal  manipulations.  It  becomes 
meaningful  when  combined  with  the  expression  for  v  and  integrated  (using  Theorem  4)  as  part  of 
the  integration  by  parts  formula.  Using  these  terms  for  u,  v,  du,  and  dv  for  integrating  the  term  in 
the  right-hand  side  of  Equation  (32)  with  the  integral, 

S — £ 

lim[  J(  0.5  -  H)(s-r)"®-^““r°^"”x(r)dr-h  e®-^'“(s-  X  (s  -  e)]  ds  = 

£-~>0  Q 


lim  [“£ 
e-^0 


0.5-H 


s-e 


^X(s-e)-f  (0.5-H)  J  (s 

0 


r)'>5-»X(r)r-''’-"dr  + 


S'”£ 

J  (s-  r)"®-^~»r  ®'^"“dX(r)  +  e®-^‘»(s-  e)®-5"”X(s  -e)]ds  = 


0 


((  0.5  -  H  )f  (s  -  r  X(r  )r  |  ^  0.5-H 

0  0 


dX(r  )]ds 


(35) 


R(t)  = - 

r(  1.5-H) 


JsH-0.5 

0 


s 

[(  0.5  -H)j  (s  -  r  X(r)r  "®-5~“dr  + 

0 


J(s-r)  “  r  ®-^-”dX(r  )]ds 
0 


(36) 
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Hence,  taking  the  derivative  of  R(t)  with  respect  to  t  gives  the  following: 


dR(t) 

dt 


^  t“-0-^[(  0.5  -  H)J(t-r)°-5-”X(r)r-®-^-»dr  + 


r(  1.5-H  ) 


J(t-r)-0-5-Hr«'5-«dX(r)] 

0 


(37) 


Following  the  same  arguments  in  Reference  2,  pp.  71-76,  for  the  derivation  for  the 
estimator  in  Equation  (19),  the  estimator  for  the  present  model  with  t  e  [0,T]  can  be  shown  to  be 
as  follows: 


T 

J(dR  (s)/ds)dY(s) 

0(T)  =  Aj - - 

J[dR(s)/ds]^ds 

0 

Values  for  Y  can  be  derived  from  Equation  (22).  Furthennore,  the  derivative  of  R  as  shown  in 
Equation  (37)  is  in  terms  of  X.  Hence,  everything  in  the  estimator  of  0  can  be  obtained  from  X, 
which  is  the  original  given  data.  This  means  that  estimator  can  be  implemented  to  give  a 
numerical  estimation  of  0. 

Now  consider  an  infinite  class  of  closed  and  bounded  sets  {[0,Tn]:  Tn  finite,  n  =  0,1,2,..., 
and  Tn+l>Tn}.  Note  that  Theorem  3  implies  that  a  Brownian  motion  exists  for  each  [0,Tn]  such 
that  the  desirable  properties  are  fulfilled.  However,  examine  the  formula  in  Theorem  3,  which 
tells  about  the  existence  of  Brownian  motion  given  FBM  in  a  finite  interval.  This  formula  is 
valid  for  both  [0,Tn]  and  [0,Tn+k]  for  any  positive  integer  k,  and  it  shows  that  the  sample  path 
up  to  Tn-hk  of  the  interval  when  restricted  to  Tn  is  the  same  as  the  sample  path  of  Brownian 
motion  derived  for  the  closed  and  bounded  set  [0,Tn].  Therefore,  the  estimator  derived  in 
Equation  (38)  is  equivalent  to  an  estimator  for  all  finite  T  e  [0,  ■»).  Thus,  considering  whether 
asymptotic  properties  such  as  consistency  exist  is  still  meaningful  in  spite  of  the  fact  that  the 
algorithm  was  derived  for  a  finite  index  set.  This  algorithm  can  be  shown  to  be  strongly 
consistent;  that  is,  it  converges  a.  s.  to  the  actual  0  by  the  same  argument  given  in  Reference  2, 
pp.  77-86,  for  showing  the  strong  consistency  of  the  parametric  estimator  for  the  model  given  by 
Equation  (19). 
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SUMMARY 


The  purpose  of  this  research  project  was  to  find  methods  for  modeling  physical 
phenomena  that  contain  an  unknown  constant  parameter  and  random  noises,  characterized  by 
long-tei-m  slowly  decreasing  time  dependencies.  The  model  being  considered  is  in  the  form  of  a 
stochastic  differential  equation  written  as  Equations  (9)  or  (10)  in  this  report.  The  unknown 
parameter  is  0.  With  this  purpose  in  mind,  an  estimator  that  estimates  the  unknown  parameter 
was  tested.  The  noise  process  being  used  in  this  project  is  fractional  Brownian  motion  (FBM  or 
Bh)  or  fractional  Gaussian  noise  (FGN  or  Wh).  The  question  concerning  whether  this  estimator 
is  convergent  (converges  to  the  actual  value  of  the  pai'ameter)  remains  inconclusive  from  the 
tests.  FBM  and  FGN  aa’e  themselves  characterized  by  a  single  parameter  known  as  H.  An 
algorithm  that  can  be  used  to  estimate  this  single  parameter  of  FBM  was  also  described  in  this 
report. 


A  Monte  Carlo  simulation  of  the  FBM  as  a  multivariate  normally  distributed  random 
vector  (where  each  component  of  the  vector  represents  FBM  at  a  different  time)  has  been 
developed.  This  was  in  turn  used  to  create  a  simulation  of  the  variable  called  X,  in  the  stochastic 
differential  equation  model,  Equation  (10).  The  simulation  of  X  was  necessary  to  test  the  results 
of  the  estimator  for  the  parameter  0.  Furthermore,  the  simulator  of  X  itself  has  potentials  for 
applications  for  studying  a  physical  process  that  cannot  be  replicated  in  a  laboratory,  although  its 
behavior  needs  to  be  studied. 

In  addition  to  the  above,  an  alternate  estimator  was  also  derived  to  estimate  the  parameter 
0.  This  estimator  must  assume  prior  knowledge  of  the  value  of  the  FBM  parameter  H  but  can  be 
shown  to  converge  to  the  true  value  of  the  parameter. 
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APPENDIX  A 

FORTRAN  LISTING  OF  THE  SIMULATOR 
AND  ESTIMATOR  PROGRAMS 


A-l/A-2 


oooo  n  n  no  n  noon  non 


I 
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PROGRAM  ESTIMATiai 

C  THIS  PROGRAM  IS  IN  DOUBLE  PRECISION  AND  IS  CALLED  ESTIMATI0N2DP.  A  SINGLE 
C  PRECISION  VERSICW  IS  IDENTICAL  WITH  THIS  PROGRAM  EXCEPT  FOR  TOE  FACT  THAT 
C  "IMPLICIT  REAL*8  (A-H,0-Z)"  IS  NOT  USED  AND  SINGLE  PRECISION  VERSICKS  OF 
C  THE  NSWC  LIBRARY  OF  MATO  ROUNTINES  ARE  USED  INSTEAD  OF  THE  DOUBLE  PRECISION 
C  VERSIONS  (REF.  8). 

C  USE  LINK  ESTIMATICN2DP,DNRVG,FILE1  FOR  LINKING  CW  THE  VAX  COMPUTER. 

IMPLICIT  REAL*8  {A-H,0-Z) 

COMMON  X(1000),BH(1000,1),U(1000) 

DIMENSICW  A( 1000, 1000), B(500500) 

C  FOR  X(IOOO),  USE  B(500500),  FOR  X(2000),  USE  8(2001000) 

PI=3. 1415926535897932 

SEE  NRVG  PROGRAM  IN  NSWC  LIBRARY  OF  MATH  ROUTINES  (REF.  8)  AND  ALIGN  THE  NEXT 
FOUR  VARIABLES  TO  CHECK  THEIR  MEANINGS.  NOTE  THAT  THIS  PROGRAM  USES  DRVG, 

THE  DOUBLE  PRECISICW  VERSIOI  OF  NRVG. 

M0=0 
NVEC=1 
KX=2000 
ISEED=2795601 

H  IS  THE  FRACTIONAL  BROWNIAN  MOTICW  (FBM)  PARAMETER,  TOETAO  IS  THE  INPUT 
PARAMETER  OF  TOE  MODEL  DX(T)-(THETA0)X(T)DT+DBH(T)  WHERE  BH  IS  FBM 
PARAMETERIZED  BY  H.  THIS  PROGRAM  SIMULATES  X(T)  AFTER  SIMULATING  BH 
VIA  THE  NRVG  ROUTINE. 

WRITE(6,*)'ENTER  H,N=#  OF  IICREMENTS,DEL=INCREMENT,TOETAO<=THETA' 
WRITE(7,*)'ENTER  H,N=#  OF  INCREMENTS, DEL=INCREMENT,THETAO=THETA' 

READ( 5,*)H,N,DEL,THETA0 

WRITE(6,*) 'H,N,DEL,THE1A0' ,H,N,DEL,THETA0 
WRITE(7,*) 'H,N,DEL,THETA0' ,H,N,DEL,THETA0 
M=N 

U  IS  THE  MEAN  VECTOR. 

DO  1=1, M 
U(I)=0. 

ENDDO 
Z=2.-2.*H 

GA»=<5AMMA  FUNCTION  (SEE  REF.  8  FOR  EXPLANATION  OF  THE  DGAMMA  FUNCTION  GIVEN 
ON  THE  NEXT  LINE). 

GAM-DGAMMA(Z) 

PIH=PI*H 
IF(H.EQ..5)raEN 
VH=1.0 
ELSE 

VH  IS  THE  SCALE  FACTOR  IN  THE  COVARIANCE  MATRIX  OF  BH. 

VH=- ( GAM*COS ( PIH ) )/( PIH* ( 2 . *H-1 . ) ) 

ENDIF 

WRITE(6,*) 'VH' ,VH 
WRITE(7,*)'VH',VH 

CALCULATE  TOE  COVARIANCE  MATRIX  OF  BH. 

DO  1=1, N 
DO  J=I,N 

TIME1=FLQAT ( I ) *DEL 
TIME2=FLOAT( J ) *DEL 
S2H=ABS(TIME1) 

S2H=S2H**(2.*H) 

T2H=ABS(TIME2) 

T2H=T2H**(2.*H) 

TS=ABS ( TIME2-TIME1 ) 

TS=TS**(2.*H) 

A( I , J ) = ( S2H+T2H-TS ) *VH/2 . 

A(J,I)=A(I,J) 

ENDDO 
ENDDO 

WRITE(6,*)  'A' 

WRITE(7,*)  'A' 

DO  1=1, N 

WRITE(6,*)  (A(I,J),  J=1,N) 


A-3 
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WRITE(7,*)  (A(I,J),  J=1,N) 

ENDDO 

DMCVFS  PACKS  THE  COVARIANCE  MATRIX  (REF  8).  THIS  IS  REQUIRED  BY  DNRVG. 

CALL  DMCVFS(A,M,M,B) 

WRITE(6,*)  'B',B 
WRITE(7,*)  'B',B 

K4RVG  SIMULATES  FBM  KNOWN  AS  BH  IN  THIS  PROGRAM  (REF  8). 

CALL  DNRVG(MO,ISEED,NVEC,M,B,U,BH,KX,IERR) 

WRITE(6,*)  'MO,ISEED,NVEC' ,MO,ISEED,NVEC 
WRITE(6,*)  'M,KX',M,KX 
WRITE(6,*) 'U' ,U 
WRITE(6,*)  'B',B 
WRITE(6,*)  'BH  VALUES',  BH 
WRITE(6,*)  'lERR  -  ',IERR 
WRITE(6,*)  'U(1),U(100)',U(1),U(100) 

WRITE(7,*)  'MO,ISEED,NVEC' ,M0,ISEED,NVEC 
WRITE(7,*)  'M,KX',M,KX 
WRITE(7,*)'U',U 
WRITE(7,*)  'B',B 
WRITE(7,*)  'BH  VALUES',  BH 
WRITE(7,*)  'lERR  =  ',IERR 
AV=0. 

DO  1=1, N 

AV=AV+BH(I,1) 

ENDDO 
C  SAMPLE  MEAN  OF  BH 
AV=AV/N 

WRITE (6,*) 'SAMPLE  MEAN  OF  BH' ,AV 
WRITE(7,*) 'SAMPLE  MEAN  OF  BH' ,AV 

C  STATE  SIMULATES  X(T)  OF  THE  PARAMETRIC  STOCHASTIC  DIFFERENTIAL  EQUATICW. 

C  ESTIM  IS  THE  PARAMETRIC  ESTIMATOR  OF  THE  STOCHASTIC  DIFF.  EQTO. 

CALL  STATE(N,THETAO,DEL) 

CALL  ESTIM(N,THETA,DEL) 

C  WRITE(6,*)  'STATE  MODEL  VALUES',  X 

WRITE  (6,*)  'ACTUAL  THETA',  THETAO 
WRITE(6,*)  'ESTIMATED  THETA',  THETA 
WRITE(7,*)  'STATE  MODEL  VALUES',  X 
WRITE  (7,*)  'ACTUAL  THETA',  THETAO 
WRITE(7,*)  'ESTIMATED  THETA',  THETA 
STOP 
END 

SUBROUTINE  STATE (N, THETAO , DEL) 

IMPLICIT  REAL*8  (A-H,0-Z) 

REAL  LNSUM,LNEXP,LNX 

COMMm  X(2000),BH(2000,1),U(2000) 

SUM=0. 

AV1=0. 

ADD=0. 

ADDSQ=0. 

C  GO  TO  10  IF  THE  ACTUAL  PARAMETER  IS  NEGATIVE. 

IF(THETAO.LT.O)GO  TO  10 
DO  1=1, N 

REALI=FLOAT( I-l ) 

REALJ=FLOAT( I ) 

IF(I.EQ.l)  THEN 
DBH=BH(1,1) 

C  WRITE(7,*) 'DBH' ,DBH 

PCWER=0 
ELSE 

DBH=BH(I,1)-BH(I-1,1) 

C  WRITE(7,*)'DBH' ,DBH 

POWER=-THETAO*REALI*DEL 

ENDIF 

C  AVI  AND  ADD  USED  TO  CALCULATE  SAMPLE  MEAN/VAR  OF  DBH  FOR  DEBUGGING  THIS  PROG 
ADD=ADD+DBH 
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C  ADDSQ  AND  ADD  USED  TO  CALC  SAMPLE  VAR.  FOR  DEBUGGING. 

DBHSQ=DBH**2 

ADDSQ=ADDSQ+DBHSQ 

FUNC=EXP(PCWER) 

TERM=FUNC*I®H 

SUM=SUM+TERM 

C  WRITE (7,*) 'SUM  FOR  LOG', SUM 

C  TAKE  LOG  IN  ORDER  TO  HELP  IN  CASE  OF  LARGE  SUM. 

C  THE  FOLLOWING  TAKES  CARE  OF  FACT  THAT  ARG.  OF  LOG  CANNOT  BE  NEG. 

IF(SUM.LT.O)  THEN 
LNSUM=LOG(-SUM) 

ELSE 

LNSUM=LOG(SUM) 

ENDIF 

LNEXP=THETA0*REALJ*DEL 

LNX=LNEXP+LNSUM 

X(I)«=EXP(LNX) 

C  NEXT  LINE  TAKES  CARE  OF  CASE  WHERE  LOG  OF  -SUM  WAS  TAKEN. 

IF(SUM.LT.O)  X(I)— X(I) 

C  ABOVE  GIVES  X(T)  AT  T=I*DEL  FOR  THE  CASE  WHERE  THETAO  IS  NONNEGATIVE. 

ENDDO 

C  WRITE(7,*)  'X  VALUE', X 

AV1=ADD/N 

WRITE (6,*) 'SAMPLE  MEAN  OF  DBH' ,AV1 
WRITE(7,*) 'SAMPLE  MEAN  OF  DBH' , AVI 
SQADD=ADD**2 

VARHAIfe  ( N*ADDSQ-SQADD )/{ N*  ( N-1 ) ) 

WRITE( 6,*) 'SAMPLE  VAR  OF  IBH' ,VARHAT 
WRITE (7,*) 'SAMPLE  VAR  OF  DBH',  VARHAT 
RETURN 

10  DO  1=1, N 
SUM=0. 

DO  J=1,I 
REALI=FLOAT( I ) 

REALJ=FLQAT(J-1) 

IF(J.EQ.l)  THEN 
DBH=BH(1,1) 

C  WRITE(7,*)'DBH',DBH 

PCWER=THETAO*REALI*DEL 

ELSE 

DBH=BH ( J , 1 ) -BH( J-1 , 1 ) 

C  WRITE(7,*)'DBH',DBH 

PCWER=THETA0* ( REALI *DEL-REALJ*DEL ) 

ENDIF 

C  AVI  AND  ADD  USED  TO  CALCULATE  SAMPLE  MEAN/VAR  OF  DBH  FOR  DEBUGGING  THIS  PROG. 
ADD=ADD+DBH 

C  ADDSQ  AND  ADD  USED  TO  CALC  SAMPLE  VAR.  EXDR  DEBUGGING. 

DBHSO=DBH**2 
ADDSQ=ADDSQ+DBHSQ 
FUNC=EXP( POWER) 

TERM=FUNC*DBH 

SUM=SUM+TERM 

ENDDO 

X(I)=SUM 

C  THIS  GIVES  X  AT  T=I*DEL  FOR  THE  CASE  WHERE  THETAO  IS  NEGATIVE. 

ENDDO 

NN=N*(N+l)/2 

AV1=ADD/NN 

WRITE (6,*) 'SAMPLE  MEAN  OF  DBH' , AVI 
WRITE (7,*) 'SAMPLE  MEAN  OF  DBH' , AVI 
SQADD=ADD**2 

C  VARHAT=(NN*ADDSQ-SQADD)/(NN*(NN-1 ) ) 

C  WRITE (6,*) 'SAMPLE  VAR  OF  DBH' , VARHAT 

C  WRITE(7,*) 'SAMPLE  VAR  OF  DBH',  VARHAT 

RETURN 
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END 

SUBRCXrriNE  ESTIM(N, THETA, DEL) 

C  THIS  IS  A  GENERALIZED  VERSION  OF  THE  LEAST  SQUARES  ESTIMATOR  FOR  THE  PARAMETER 
C  THETA  IN  A  STOCHASTIC  DIFFERENTIAL  EQUATION  DX(T)=(THETA)X(T)DT+DBH(T) .  SEE 
C  THE  TEXT  OF  THIS  REPORT  FOR  AN  EXPLANATION. 

IMPLICIT  REAL*8  (A-H,0-Z) 

COMMCXa  X(2000),BH(2000,1),U(2000) 

SUMl-0. 

SUM2=0. 

DO  I«1,N 

IF(I.EQ.1)THEN 

XTERM1=0. 

DX=X(I) 

ELSE 

DX=X(I)-X(I-1) 

XTERM1=X( I-l ) *DX 
ENDIF 

SUM1=SUM1+XTERM1 

C  SUHl  IS  THE  INTEGRAL  OF  THE  PROCESS  X  WITH  RESPECT  TO  THE  PROCESS 
C  X  ITSELF.  INTEGRATION  IS  FROM  0  TO  TIME  N*DEL. 

ENDDO 
DO  J=1,N 

IF(J.EQ.1)THEN 

XTERM2=0. 

ELSE 

XSQ=X{ J-l)**2. 

XTERM2=XSQ*DEL 

ENDIF 

C  SUM2  IS  THE  INTEGRAL  OF  THE  SQUARE  OF  THE  PROCESS  X  WITH  RESPECT  TO  TIME. 

C  INTEGRATION  IS  FROM  0  TO  TIME  N*DEL. 

SUM2=SUM2+XTEBM2 
ENDDO 

WRITE(6,*)'SUM1' ,SUM1 
WRITE(6,*)'SUM2' ,SUM2 
WRITE( 7 , * ) ' SDMI ' , SUMl 
WRITE( 7 , * ) ' SUM2 ' , SUM2 
THETA=SUM1/SUM2 

NOTE  THAT  IN  BOTH  INTEGRALS  RESPRESENTED  BY  SUMl  AND  SUM2  RESPECTIVELY 
ARE  APPROXIMATED  WHERE  IN  EACH  TERM  OF  THE  SUM,  THE  INTEGRAND  IS  APPROX. 

AS  THE  I-l  TERM.  THIS  IS  TO  MIMIC  THE  ITO  INTEGRAL.  SINCE  THESE  ARE 
NOT  TRUE  ITO  INTEGRALS,  IT  IS  NOT  NECESSARY  TO  CHOOSE  THE  TERM  AT  I-l. 

RETURN 
END 
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